## Shorter R code to run models and get tables
## For the Robustness check
## Appendix: using alternative measures of threat which prioritize the likelihood


## 2. Main models -------------------------------------------------------------
## Fixed effect models
fe1 <- lm(LJI ~ threatS + country_name, dat); summary(fe1)
fe2 <- lm(LJI ~ threatS  + gdpS + demdurS + country_name,dat); summary(fe2)
## IV models
iv_total1 <- ivreg(LJI ~ threatS + gdpS + demdurS | gdpS +  demdurS + ivvar1,
                   data=dat); summary(iv_total1)  
iv_total2 <- ivreg(LJI ~ threatS + gdpS + demdurS | gdpS +  demdurS  +  ivvar1 + ivvar2,
                   data=dat); summary(iv_total2)    

## Robust standard errors for fixed effects models ---------------------------
G <- length(unique(dat$country_name))
N <- nobs(fe1)
dfa <- (G/(G-1))* (N-1)/fe1$df.residual
clustered_vcov <- dfa * vcovHC(fe1, type="HC0", cluster="group", adjust=TRUE)
G <- length(unique(dat$country_name))
N <- nobs(fe2)
dfa <- (G/(G-1))* (N-1)/fe2$df.residual
clustered_vcov2 <- dfa * vcovHC(fe2, type="HC0", cluster="group", adjust=TRUE)
fe1.robust.se <- sqrt(diag(clustered_vcov))
fe2.robust.se <- sqrt(diag(clustered_vcov2))
fe1.coeftest <- coeftest(fe1, vcov=clustered_vcov)
fe1.pvals <- fe1.coeftest[,4]
fe2.coeftest <- coeftest(fe2, vcov=clustered_vcov2)
fe2.pvals <- fe2.coeftest[,4]

## Robust se for IV models -----------------------------------------------------
# texreg does not by default include IV gof statistics
G <- length(unique(dat$country_name))
N <- nobs(iv_total1)
dfa <- (G/(G-1))* (N-1)/iv_total1$df.residual
clustered_vcov_ivtotal1 <- dfa * vcovHC(iv_total1, type="HC0", cluster="group", 
                                        adjust=TRUE)
V1.robust.se <- sqrt(diag(clustered_vcov_ivtotal1))
V1.coeftest <- coeftest(iv_total1, vcov= clustered_vcov_ivtotal1)
V1.pvals <- V1.coeftest[,4]
tr <- texreg::extract(iv_total1,include.rsquared=FALSE,include.adjrs=FALSE)
tr@gof <- c(tr@gof, summary(iv_total1,diagnostics=TRUE)$diagnostics[1,4])
tr@gof.names <- c(tr@gof.names, "Weak instruments test")
tr@gof.decimal <- c(tr@gof.decimal, TRUE)
tr@gof <- c(tr@gof, summary(iv_total1,diagnostics=TRUE)$diagnostics[2,4])
tr@gof.names <- c(tr@gof.names, "Wu-Hausman test")
tr@gof.decimal <- c(tr@gof.decimal, TRUE)
tr@gof <- c(tr@gof, summary(iv_total1,diagnostics=TRUE)$diagnostics[3,4])
tr@gof.names <- c(tr@gof.names, "Sargan test")
tr@gof.decimal <- c(tr@gof.decimal, TRUE)
## iv total 2 
N <- nobs(iv_total2)
dfa2 <- (G/(G-1))* (N-1)/iv_total2$df.residual
clustered_vcov_ivtotal2 <- dfa2 * vcovHC(iv_total2, type="HC0", cluster="group", 
                                         adjust=TRUE)
V2.robust.se <- sqrt(diag(clustered_vcov_ivtotal2))
V2.coeftest <- coeftest(iv_total2, vcov= clustered_vcov_ivtotal2)
V2.pvals <- V2.coeftest[,4]
tr2 <- texreg::extract(iv_total2,include.rsquared=FALSE,include.adjrs=FALSE)
tr2@gof <- c(tr2@gof, summary(iv_total2,diagnostics=TRUE)$diagnostics[1,4])
tr2@gof.names <- c(tr2@gof.names, "Weak instruments test")
tr2@gof.decimal <- c(tr2@gof.decimal, TRUE)
tr2@gof <- c(tr2@gof, summary(iv_total2,diagnostics=TRUE)$diagnostics[2,4])
tr2@gof.names <- c(tr2@gof.names, "Wu-Hausman test")
tr2@gof.decimal <- c(tr2@gof.decimal, TRUE)
tr2@gof <- c(tr2@gof, summary(iv_total2,diagnostics=TRUE)$diagnostics[3,4])
tr2@gof.names <- c(tr2@gof.names, "Sargan test")
tr2@gof.decimal <- c(tr2@gof.decimal, TRUE)
## save the table  ------------------------------------------------------------
texreg(list(fe1, fe2, tr, tr2),
       file = dir_table,
       override.se = list(fe1.robust.se, fe2.robust.se, V1.robust.se, V2.robust.se),
       override.pvalues= list(fe1.pvals, fe2.pvals, V1.pvals, V2.pvals),
       caption= texcaption,
       omit.coef="name",caption.above=TRUE,
       booktabs=TRUE,dcolumn=TRUE,use.packages=FALSE,
       label=texlabel,
       custom.coef.names=c("Intercept","Demand for insurance","log(GDP/capita)",
                           "Years democratic (logged)"),
       include.aic=FALSE,include.deviance=FALSE,
       include.rmse=FALSE,style="apsr",digits=2,scalebox=0.75, float.pos="t!") 

